#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <list>


#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
//#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;

///---------------------------------------------------------------------
/// This is a R-K 4 version for particle traking with respect to vortex
///
#include "ParTrack/RungeKutta4Vort.hpp"
///---------------------------------------------------------------------

void writeTofile_DX(list<floatTinyArray_t> &, int, string);
inline double dot(floatTinyArray_t &, floatTinyArray_t &);

int main () 
{
///
/// Mesh definition
///
  StructuredMesh< Uniform<double, 2> > mesh(1, 257, 1, 257);
    floatTinyArray_t deltas = mesh.getDeltas();
///
    double dt, tol, dist, tolf, fold;

    ofstream flong("longitud");

///
/// Middle line construction
///
    int Npoints, Nsteps, frec, Ncicles ;
    cout << " Npoints = "; cin >> Npoints;
    cout << " Steps = "; cin >> Nsteps;
    cout << " Time step = "; cin >> dt;
    cout << " Tolerance = "; cin >> tol;
    cout << " Fold = "; cin >> tolf;
    cout << " No Cicles = "; cin >> Ncicles;
    cout << " Print freq = "; cin >> frec;

///---------------------------------------------------------------------
/// This part is for reading the vortex trayectory of 2D chaotic mixing
///
    ifstream ftrayvort("trayZero.dat");

    list<floatTinyArray_t> trayvort;
    list<floatTinyArray_t>::iterator iter_v;
    floatTinyArray_t vort_pos, vel_vort;
    for(int i = 0; i <= 801; ++i) {
	ftrayvort >> vort_pos(0) >> vort_pos(1);
	trayvort.push_back(vort_pos);
    }


//    for(iter_v = trayvort.begin(); iter_v != trayvort.end(); ++iter_v) {
//	cout << " (x,y) = (" << (*iter_v)(0) << ", " << (*iter_v)(1) << ")\n";
//    }
///---------------------------------------------------------------------

    list<floatTinyArray_t> linec;
    list<floatTinyArray_t>::iterator iter, iter_p1;

    floatTinyArray_t aux, x_i_n, x_i_np1, x_im1_n, x_im1_np1, 
	x_i, x_im1, x_ip1;
    floatTinyArray_t x_ip1_n, r_im1, r_ip1;
    double left_extreme = deltas(0);
    double right_extreme = 1.0 - deltas(0);

    double factor = (right_extreme - left_extreme) / (Npoints - 1);

    for(int i = 0; i < Npoints; i++) {
	aux(0) = left_extreme + i * factor;
	aux(1) = 0.5;
	linec.push_back(aux);
    }

    writeTofile_DX(linec, 0, "line.");
    cout << "\n ----- Begining track ... \n\n";
///
/// Velocity field at time t
///
    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());
///
/// Velocity field at time t + dt
///
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());
///
/// Particle definition
///
    ParticleTracking< RungeKutta4< Linear <double, 2 > > > particle;
///
/// Follow all the particles for several cicles...
///
    int count = 1;
    for(int cicles = 1; cicles <= Ncicles; ++cicles) {
	cout << " ----> Cicle " << cicles << endl;
///---------------------------------------------------------------------
/// This is a R-K 4 version for particle traking with respect to vortex
///
	iter_v = trayvort.begin();
	cout << "primer valor = " << (*iter_v) << endl;
///---------------------------------------------------------------------
	for(int step = 1; step <= Nsteps; ++step) {
///
/// Read the velocity at time t from a file
///
	    InOut::readFromFile_1(u, v, step, "velc.");
///
/// Read the velocity at time t + dt from a file
///
	    InOut::readFromFile_1(un, vn, step + 1, "velc.");
///
/// Track every particle on the line
///
/// Initial point in the line  x(i-1)^n
	    iter = linec.begin(); 
	    x_im1_n = *iter;      
/// Track first point to the next position x(i-1)^(n+1) 
/////	    *iter = particle.track(*iter, deltas, u, v, un, vn, dt);

///---------------------------------------------------------------------
/// This is a R-K 4 version for particle traking with respect to vortex
///
	    vel_vort = *iter_v;
	    iter_v++;

	    cout << " ---> x_1 = " << vel_vort << endl;
	    cout << " ---> x_2 = " << (*iter_v) << endl;
	    vel_vort(0) = ( (*iter_v)(0) - vel_vort(0) ) / dt;
	    vel_vort(1) = ( (*iter_v)(1) - vel_vort(1) ) / dt;

	    cout << " ---> vel_vort = " << vel_vort << endl;
	    
	    *iter = particle.track(*iter, deltas, u, v, un, vn, dt, vel_vort);
///---------------------------------------------------------------------

	    x_im1_np1 = *iter; 
/// Advance iter to the second point in the line
	    iter++;  

	    for( ; iter != linec.end(); ++iter) {
/// Current point  x(i)^n
		x_i_n = *iter; 
/// Track the current point to the next position, x(i)^(n+1)
/////		x_i_np1 = particle.track(*iter, deltas, u, v, un, vn, dt);

///---------------------------------------------------------------------
/// This is a R-K 4 version for particle traking with respect to vortex
///
		x_i_np1 = particle.track(*iter,deltas,u,v,un,vn,dt,vel_vort);
///---------------------------------------------------------------------


/***

/// Calculate the distance between consecutive tracked points 
/// x(i-1)^(n+1) and x(i)^(n+1)
		dist = sqrt( ( x_i_np1(0) - x_im1_np1(0) ) * 
			     ( x_i_np1(0) - x_im1_np1(0) ) +
			     ( x_i_np1(1) - x_im1_np1(1) ) * 
			     ( x_i_np1(1) - x_im1_np1(1) ) );
/// if dist exceeds the tolerance, adds a point between points: 
/// x(i)^n and x(i-1)^n
		
		if ( dist > tol ) {
		    aux(0) = ( x_i_n(0) + x_im1_n(0) ) * 0.5;
		    aux(1) = ( x_i_n(1) + x_im1_n(1) ) * 0.5;
		    linec.insert(iter,aux);
/// Rewind the iterator 2 positions, because the for will advance one
/// position.
		    iter--; iter--;
		} else {
/// If dist does not exceeds the tolerance, then re-state the value
/// of x(i-1)^n, x(i-1)^(n+1) and modify the new position of x(i) in the
/// list ( x(i)^(n+1) )
		    x_im1_n   = x_i_n;
		    x_im1_np1 = x_i_np1;
		    *iter     = x_i_np1;
		}
***/
	    }

	    if( !(count % frec) ) {


/***
		iter = linec.begin();
		x_im1 = *iter;
		iter++;
		for( ; iter != --linec.end(); ++iter) {
		    x_i = *iter;
		    iter_p1 = iter;
		    x_ip1 = *(++iter_p1);
		    r_ip1(0) = x_i(0) - x_ip1(0);
		    r_ip1(1) = x_i(1) - x_ip1(1);
		    r_im1(0) = x_i(0) - x_im1(0);
		    r_im1(1) = x_i(1) - x_im1(1);
		    fold = 0.5 * dot(r_im1, r_ip1) +  0.5;
		    if ( fold <= tolf ) {
			iter_p1 = iter;
			iter--;
			linec.erase(iter_p1);
		    } else {
			x_im1 = x_i;
		    }
		}
***/
		writeTofile_DX(linec, count, "line.");


	    }

	    flong << count << "\t" << linec.size() << endl;
	    
	    cout << " step = " << step 
		 << " count = " << count 
		 << " Npoints = " << linec.size()
		 << endl;


	    count++;
	}


    }
    
    cout << "\n End ... \n\n";

    flong.close();

    return 0;
}

inline double dot(floatTinyArray_t &r_im1, floatTinyArray_t &r_ip1) 
{
    double mm1 = sqrt( r_im1(0) * r_im1(0) + r_im1(1) * r_im1(1) );
    double mp1 = sqrt( r_ip1(0) * r_ip1(0) + r_ip1(1) * r_ip1(1) );
    return ( (r_im1(0) * r_ip1(0) + r_im1(1) * r_ip1(1)) / ( mm1 * mp1 ) );
}


void writeTofile_DX(list<floatTinyArray_t> &linec, int n, string filename)
{

    int Npoints = linec.size();
    list<floatTinyArray_t>::iterator iter;

    ostringstream extnumber;
    extnumber.width(5);
    extnumber.fill('0');
    extnumber << n;
    string snumber = extnumber.str();
    filename += snumber;
    ofstream file (filename.c_str());  

    cout << "writing " << filename << endl;
    
    file << " object 1 class array type float rank 1 shape 2 items "
	 << Npoints << " data follows\n";
    for(iter = linec.begin(); iter != linec.end(); ++iter)
	file << (*iter)(0) << "   " << (*iter)(1) << "   ";
    file << endl;

    file << "object 2 class array type int rank 1 shape 2 items "
	 << Npoints - 1 << " data follows\n";
    for(int i = 0; i < Npoints - 1; ++i)
	file << i << "   " << i + 1 << "   ";
    file << endl;

    file << "attribute \"element type\" string \"lines\"\n"
	 << "attribute \"ref\" string \"positions\"\n";
    file << "object 3 class array type float rank 0 items "
	 << Npoints << " data follows\n";
    for(int i = 0; i < Npoints; ++i)
	file << 1 << "\t";
    file << endl;
 
    file << "attribute \"dep\" string \"positions\"\n";
    file << "object \"irregular positions regular connections\" class field\n"
	 << "component \"positions\" value 1\n"
	 << "component \"connections\" value 2\n" 
	 << "component \"data\" value 3\n" 
	 << "end\n" ;
    file.close();

}


